/*! 
 ***************************************************************************
 * \file singlePhase3D.cpp
 * Single phase in porous media in 1D.
 *
 * \par Description
 * This code solves the 3D single phase equation for pressure:
 * \f$ \frac{\phi \mu c}{k}\frac{\partial p}{\partial t} = \nabla^2 p \f$.
 * This equation is solved in \f$ (x,y,z) \in [0,100 cm] \times [0,10 cm] \times [0,10 cm] \f$.
 * The data used are: 
 * The initial condition for the temperature is \f$ P = 0 \f$
 * and the boundary conditions:  \f$ P_L=2 atm \f$ and \f$ P_R=1 atm\f$.
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./singlePhase  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ Mon Mar 17 08:47:37 GMT 2008 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,3>::huge ScalarField3D;


void exactSolution(ScalarField3D &T, double hx, double dx, double PL, double PR)
{
  int NI = T.ubound(firstDim);
  int NJ = T.ubound(secondDim);
  int NK = T.ubound(thirdDim);
  
  for(int i = 1; i < NI; ++i)
    for(int j = 1; j < NJ; ++j) 
      for(int k = 1; k < NK; ++k)
	T(i,j,k) = PL + (PR - PL) * i * dx / hx;
}

int main() 
{
  double length_x, length_y, length_z, dt, tolerance = 1e-7; 
  int num_nodes_x, num_nodes_y, num_nodes_z, frequency, tdma_iter = 5000;
  double porosity, viscosity, permeability, compresibility;
  double PL, PR, Tmax;

// Reading data from the input file
    std::ifstream input_cin ("inputSingle3D"); 
    input_cin >> length_x >> length_y >> length_z
	      >> num_nodes_x >> num_nodes_y >> num_nodes_z
	      >> porosity
	      >> permeability
	      >> viscosity
	      >> compresibility
	      >> PL
	      >> PR
	      >> dt
	      >> Tmax
	      >> frequency;
    input_cin.close();

    Tmax += dt;

    double Gamma = permeability / (porosity * viscosity * compresibility);

    std::cout << "\n Gamma = " << Gamma;

    StructuredMesh<Uniform<double, 3> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y,
					     length_z, num_nodes_z);
    double dx = mesh.getDelta(X), dy = mesh.getDelta(Y), dz = mesh.getDelta(Z);
    mesh.print();
    ScalarField3D p ( mesh.getExtentVolumes() );
    ScalarField3D pa( mesh.getExtentVolumes() );
    ScalarField3D pn( mesh.getExtentNodes()   );

    p = 1.0; // initial condition on volumes
    pa = 1.0; // initial condition on nodes
    pn = 1.0; // initial condition on nodes

    Range all = Range::all();
    p(0, all, all) = 2.0;
    pa(0, all, all) = 2.0; // boundary condition
    pn(0, all, all) = 2.0; // boundary condition

    exactSolution(pn, length_x, dx, PL, PR);
    NumUtils::interpolateToVolumes(pa, pn);

    InOut::writeToFile_DX(pn, 0, "./DataSingle3D/pres.", dx, dy, dz);
       
    SparseMatrix< Diagonal< double, 3> > A(num_nodes_x, num_nodes_y, num_nodes_z); 
    ScalarField3D              b(num_nodes_x, num_nodes_y, num_nodes_z);

    ScalarEquation< CDS<double, 3> > single_phase(p, A, b, mesh.getDeltas());
    single_phase.setDeltaTime(dt); 
    single_phase.setGamma(Gamma);
    single_phase.setDirichlet(LEFT_WALL); 
    single_phase.setDirichlet(RIGHT_WALL);
    single_phase.setNeumann(TOP_WALL);
    single_phase.setNeumann(BOTTOM_WALL);
    single_phase.setNeumann(FRONT_WALL);
    single_phase.setNeumann(BACK_WALL);
    single_phase.print();

    ofstream error_fp("./DataSingle3D/error"), residual_fp("./DataSingle3D/residual");
    int iteration = 0;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;

    while (t <= Tmax) {
      std::cout << "\n ---- Iteration " << ++iteration 
		<< "\t time = " << t
		<< "\t Res = " << residual
		<< "\t Err = " << error;
      single_phase.calcCoefficients();
      Solver::TDMA3D(single_phase, tolerance, tdma_iter, 1.0);
      residual = single_phase.calcResidual();
      error = NumUtils::calcErrorL1( p, pa );    
      //single_phase.calcErrorL1();

      error_fp << iteration << "\t" << error << std::endl;
      residual_fp << iteration << "\t" << residual << std::endl;	
 
      single_phase.update();

      t += dt;

      if ( !(iteration % frequency) ) {
	NumUtils::interpolateToNodes(pn, p);
	InOut::writeToFile_DX(pn, iteration, "./DataSingle3D/pres.",dx, dy, dz);	
      }
    }

    std::cout << "\n\n NOTE: The data were written in the [./DataSingle3D] directory \n\n";

    return 0;
}




